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Abstract: We investigate variants of Goddard's problems for nonvertical trajectories. The 
control is the thrust force, and the objective is to maximize a certain final cost, typically, 
the final mass. In this report, performing an analysis based on the Pontryagin Maximum 
Principle, we prove that optimal trajectories may involve singular arcs (along which the 
norm of the thrust is neither zero nor maximal), that are computed and characterized. 
Numerical simulations are carried out, both with direct and indirect methods, demonstrating 
the relevance of taking into account singular arcs in the control strategy. The indirect 
method we use is based on our previous theoretical analysis and consists in combining a 
shooting method with an homotopic method. The homotopic approach leads to a quadratic 
regularization of the problem and is a way to tackle with the problem of nonsmoothness of 
the optimal control. 

Key-words: Optimal control, Goddard's problem, singular trajectories, shooting method, 
homotopy, direct methods. 



Support from the French space agency CNES (Centre National d'Etudes Spatial) is gratefully acknowl- 
edged. 

* CMAP UMR CNRS 7641, Ecole Polytechnique, and INPJA Futurs, 91128 Palaiseau, France 
T Frederic . BonnansOpolytechnique . edu 
^ martinon@cmap.polytechnique.fr 

§ Universite d'Orleans, Math., Labo. MAPMO, UMR CNRS 6628, Route de Chartres, BP 6759, 45067 
Orleans cedex 2, France (Emmanuel.Trelatauniv-orleans.fr) 



Unite de recherche INRIA Futurs 
Pare Club Orsay Universite, ZAC des Vignes, 
4, rue Jacques Monod, 91893 ORSAY Cedex (France) 

Telephone : +33 1 72 92 59 00 — Telecopie : +33 1 60 19 66 08 



Arcs singuliers dans le probleme de Goddard generalise 



Resume : On etudie des variantes du probleme de Goddard pour des trajectoires non 
verticales. Le controle est la force de poussee, et l'objectif est de maximiser un certain 
critere, typiqucmcnt la masse finale. Dans ce rapport, on prouve par une analyse basee sur 
le Principe du Maximum de Pontryagin que les trajectoires optimales peuvent comporter 
des arcs singuliers (sur lesquels le module de la poussee n'est ni mil ni maximal), que Ton 
peut caracteriser et calculer. On presente des simulations numeriques, a la fois avec des 
methodes directes et indirectes, qui montrent la pertinence de la prise en compte des arcs 
singuliers dans la strategie de controle. La methode indirecte utilisee est basee sur F analyse 
theorique mentionnce preccdcmmcnt, et consiste a combiner une methode de tir avec une 
approche homotopique. L'homotopie prend la forme d'une regularisation quadratique du 
probleme, et permet de traiter l'aspect non lisse du controle optimal. 

Mots-cles : Controle optimal, probleme de Goddard, trajectoires singulicres, methode de 
tir, homotopic, methodes directes. 



Example of RR.sty 



3 



1 Introduction 

The classical Goddard's problem (see [THl [2D1 122] ) consists in maximizing the final altitude 
of a rocket with vertical trajectory, the controls being the norm and direction of the thrust 
force. Due to nonlinear effects of aerodynamic forces, the optimal strategy may involve 
subarcs along which the thrust is neither zero nor equal to its maximal value, namely, since 
the control variable enters linearly in the dynamics and the cost function is over the final 
cost, singular arcs. A natural extension of this model for nonvertical trajectories is the 
control system 

f = v, 

D(r,v) v u 

v = Tl\~9{r) + C— , 1 

m \\v\\ m 

m = — b\\u\\, 

where the state variables are r(t) £ R 3 (position of the spacecraft), v(t) £ R 3 (velocity 
vector) and m(t) (mass of the engine). Also, D(r, v) > is the drag component, g(r) £ M 3 
is the usual gravity force, and & is a positive real number depending on the engine. The 
thrust force is Cu(t), where C > is the maximal thrust, and the control is the normalized 
thrust u(t) £ IR 3 , submitted to the constraint 

||u(f)|| < 1. (2) 

The real number b > is such that the speed of ejection is C/b. Here, and throughout the 
paper, || || denotes the usual Euclidean norm in IR 3 . 

We consider the optimal control problem of steering the system from a given initial point 

r(0) = r , v(0) = v , m(0) = m , (3) 

to a certain target Mi C R 7 , in time tf that may be fixed or not, while maximizing a final 
cost. For the moment, there is no need to be more specific with final conditions and the 
cost. In real applications, the problem is typically to reach a given orbit, either in minimal 
time with a constraint on the final mass, or by maximizing the final mass, or a compromise 
between the final mass and time to reach the orbit. In our numerical experiments we will 
study the problem of maximizing the final mass (i.e., minimizing the fuel consumption) 
subject to a fixed final position r(tf) = rf, the final velocity vector and final time being 
free. 

Depending on the features of the problem (initial and final conditions, mass/thrust ratio, 
etc), it is known that control strategies that consist in choosing the control so that ||it(t)|| 
is piecewise constant all along the flight, either equal to or to the maximal authorized 
value 1, may not be optimal, as a consequence of the high values of the drag for high speed. 
Optimal trajectories may indeed involve singular arcs, and it is precisely the aim to this 
article to perform such an analysis and prove that the use of singular arcs is relevant in the 
problem of launchers. 



RR n° 0123456789 



4 



Bonnans, Martinon & Trelat 



The article is structured as follows. In Section [21 we recall the Pontryagin Maximum 
Principle, and the concept of singular trajectories. A precise analysis of the optimal control 
problem is performed in Section[3l where extremals are derived, and singular trajectories are 
computed. Theorem [T] makes precise the structure of the optimal trajectories. Section [?] is 
devoted to numerical simulations. The problem is first implemented with indirect methods, 
based on our theoretical analysis with the maximum principle, and, numerically, our method 
uses a shooting method combined with an homotopic approach. The homotopic method, 
leading to a quadratic regularization, permits to tackle with the problem of nonsmoothness 
of the optimal control. Experiments are also made using direct methods, i.e., by discretizing 
control variables and solving the resulting nonlinear optimization problem. Less precise than 
the indirect one, this method permits however to validate our approach by checking that 
results are consistent with the previously computed solution. 

Our results show, as expected, that taking into account singular arcs in the control 
strategy permits to improve slightly the optimization criterion. The numerical simulations 
presented in this paper, using a simplified and more academic model and set of parameters, 
constitute the first step in the study of a realistic launcher problem. 

2 Preliminaries 

In this section we recall a general version of the Pontryagin Maximum Principle (see [16j . 
and for instance 7 for its practical application), and a definition and characterizations of 
singular arcs. 

Consider the autonomous control system in R" 

i(t) = /(z(t), «(«)), (4) 

where / : 1R x H™ x M m — > H™ is of class C , and where the controls are measurable and 
bounded functions defined on a subintcrval [0, t e {u){ of R + with values in C H m . Let Mq 
and Mi be subsets of H™. Denote by U the set of admissible controls u, whose associated 
trajectories are well defined and join an initial point in Mq to a final point in M±, in time 
t(u) < t e (u). 

Define the cost of a control u on [0, t] by 

C(t, u) = [ f(x(s),u(s))ds + g°(t, *(*)), 
Jo 

where f° : E" x R m — ► H and g° : R" -> 1R are of class C 1 , and x(-) is the trajectory 
solution of (JTJ) associated to the control u. 

Consider the optimal control problem of finding a trajectory joining Mq to M\ and 
minimizing the cost. The final time may be free or not. 
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2.1 Pontryagin Maximum Principle 

According to the Pontryagin Maximum Principle (see [16]), if the control u £lA associated to 
the trajectory x(-) is optimal on [0, T], then there exists an absolutely continuous mapping 
p(-) : [0,1*] — ► R™ called adjoint vector, and a real number p° < 0, such that the couple 
(p(-),p°) is nontrivial, and such that, for almost every t S [0, T], 

OH 

x(t) = — (x(t),p(t),p°,u(t)), 

dH (5) 
p(t) = - — (x(t),p{t),p°,u(t)), 

where H(x,p,p° ,u) = (p, f(x,u)) + p°f°(x,u) is the Hamiltonian of the optimal control 
problem. Moreover, the function 

1 1 — > max H(x(t),p(t),p , v) 

is constant on [0,T], and the maximization condition 

H{x(t),p(t),p°,u{t)) = m & xH{x(t),p(t),p°,v) (6) 

holds almost everywhere on [0, T\. 

Moreover, if the final time T to join the target set Mi is free, then 

m^H(x(t),p(t),p°,v) = - p o?£-(T t x(T)). (7) 

for every t £ [0,T]. 

Furthermore, if Mq and Mi (or just one of them) are submanifolds of M" having tangent 
spaces in x(0) € Mq and x(T) £ Mi, then the adjoint vector can be chosen so as to satisfy 
the transversality conditions at both extremities (or just one of them) 

p(0) 1 T x(0) Mo (8) 

and 

p(T)-p°^-(T,x(T)) 1 T x{T) Mi. (9) 

An extremal of the optimal control problem is a fourth-tuple (x(-),p(-),p° ,u(-)) solution of 
([!]) and If po = 0, then the extremal is said to be abnormal, and if p° ^ then the 
extremal is said to be normal. 
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2.2 Singular arcs 

Given xq £ H" and two real numbers to, t\, with to < t\, denote by U XOt t ,ti the set of 
controls u S L°°([io, fi], 01), with 01 an open subset of O, such that the trajectory t \— * 
x(t, xq, to, u), solution of ([T]), associated with the control u on [to, fi], and such that x(to) = 
Xo, is well defined on [to,ti]- Define the end-point mapping E Xo _ tat t t by E XOtta .t 1 {u) '■— 
x(ti,Xo,to,u), for every u e U XqM m- ^ is classical that E XqMM : U XQ<to M ^ H™ is a 
smooth map. 

A control it € Ux 0i t ,ti is said to be singular if u is a critical point of the end-point 
mapping E XQ ^ ,t 1 , i.e., its differential dE Xo j ,tA u ) a ^ u is not surjective. In this case, the 
trajectory x{-,xo,to,u) is said to be singular on [to,ti]- 

Recall the two following standard characterizations of singular controls (see 0HE]). A 
control u 6 Ux ,t a ,u is singular if and only if the linearized system along the trajectory 
x(-, xq, to, u) on [tojti] is not controllable. This is also equivalent to the existence of an 
absolutely continuous mapping p\ : [to,ti] — > R n \{0} such that, for almost every t £ [to,ti], 

m = -^-(x(t), Pl {t), U (t)), m = —^±(t, x (t), P1 (t),u(t)), 

OH 

—±(x(t), Pl (t),u(t)) = 0, 
ou 

where H\(x,pi, u) = (pi, f(x, u)) is the Hamiltonian of the system. 

Note that singular trajectories coincide with projections of abnormal extremals for which 
the maximization condition ([6]) reduces to jS- = Q. 

For a given trajectory x(-) of the system ^ on [0, T], associated to a control u S M x (o),o.T: 
we say that x(-) involves a singular arc, defined on the subinterval [to,ti] C [0, T], whenever 
the control U|[ to tl ] for the control system restricted to [to,ti] is singular. 

In the case when the dynamics / and the instantaneous cost / are linear in the control u, 
a singular arc corresponds to an arc along which one is unable to compute the control directly 
from the maximization condition of the Pontryagin maximum principle (at the contrary of 
the bang-bang situation). Indeed, in this case, the above condition = along the 
arc means that some function (called switching function) vanishes identically along the arc. 
Then, it is well known that, in order to derive an expression of the control along such an arc, 
one has to differentiate this relation until the control appears explicitly. It is also well known 
that such singular arcs, whenever they occur, may be optimal. Their optimal status may be 
proved using generalized Legendre-Clebsch type conditions or the theory of conjugate points 
(see [IB [U], or see for a complete second-order optimality theory of singular arcs). 
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3 Analysis of the optimal control problem 

With respect to the notations used in the previous section, we set 

v 

x = ( v ] e R 3 x R 3 x R, f(x,u) = ( -^^| - g(r) + C-; 

-b\\u\\ 

and f° — 0. Here, the set fl of constraints on the control is the closed unit ball of 1R 3 , 
centered at 0. 

Consider the optimal control problem of minimizing some final cost g°(tf,x(tf)), for the 
control system (TTJ), with initial conditions © and final conditions x(tr) £ Mi in time tf 
which may be free or not. 

We make the following assumption. 

Assumption (H). The function g° is such that: 

• either the final mass m(tf) is free, and ^ 0, 

• or the final time tf is free, and ^ 0. 

In the first situation, the target set M\ C R 7 can be written as M\ = N% X K, where 
Ni is a subset of R 6 . A typical example is the problem of maximizing the final mass, for 
which g°(t,x) = —m. If the final condition is r(i/) = r\ and ||u(£/)|| = a, then M\ = 
{ri} x 5(0, a) x R, where 5(0, a) is the sphere of R 3 , centered at 0, with radius a. 

In the second situation, a typical example is the minimal time problem to reach some 
target. In this case, g°(t, x) = t. 

3.1 Computation of extremals 

According to Section 12. 1[ the Hamiltonian of the optimal control problem under considera- 
tion is 

H = ( Pr ,v) + (p v , _Eh^JL _ g{r) + C ±) - Pm b\\u\\, (10) 
\ m \\v\\ m I 

where ( , ) denotes the usual scalar product in R 3 . Here, the adjoint vector is denoted by 

(Vr{t)\ 

P(t) = Pv(t) G R 3 x R 3 x R. 

\Pm(t) 
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In what follows, we assume the mappings D and g to be of class C . Applying Pontryagin's 
Maximum Principle leads to the adjoint equations 

. = 1 (pv,v) 3D / dg 
r m \\v\\ dr \ ' dr 

1 (pv,v) 3D D p v D v 

Pv = ~Pr + fj Tj 7J I TT - !? (P«' U )|ni3' ( U > 

m \\v\\ ov m \\v\\ m \\ v \\ 
1 / D(r,v) v _u 

Prn = — (Pv, Tj-jf +6 — 

m \ m \\v\\ m 

Moreover, if u is an optimal control on [0, tf] , then, for almost every t £ [0, tf] , u(t) maximizes 
the function 

C 

~ —7-r{pv(t),w) - bp m (t)\\w\\, 

m{t) 

among all possible w 6 R 3 such that ||w|| < 1. 

The next technical lemma is the first step in the analysis of extremals. 

Lemma 3.1. If there exists to £ [0, tf] such that p r (to) = p v {to) = 0, then p r (t) = p v (t) = 0, 
and p m (t) = p m (tf), for every t £ [0,t/]. Moreover, p m (tf) 7^ 0, and if p m (tf) > then 
u(t) = on [0,tf], otherwise \\u(t)\\ = 1 on [0,t/]. 

Proof. The first statement follows immediately from a uniqueness argument applied to the 
system (TTTT) . It follows from the expression of the Hamiltonian function that, if p m (t) > 0, 
then u(t) = 0, and if p m (t) < 0, then ||u(f)|| = 1. In the first case of Assumption (H), the 
transversality condition ^ yields in particular 

P m (tf)=P°^(t f ,x(t f )). 

Therefore, p m (t) cannot be equal to zero (otherwise the adjoint vector (p,p°) would be zero, 
contradicting the maximum principle). In the second case of Assumption (H), it follows 
from ||7J) and (JTUJ) that 

da 

Pm{t)b\\u{t)\\=p°-^{tf,x{t } )). 
Therefore, similarly, p m {t) cannot be equal to zero. The conclusion follows. □ 



An extremal satisfying the conditions of Lemma 13.11 (ie p r (t) = p v (t) = for every 
t £ [0,£/]) is called degenerate. For such extremals, the control is either identically equal 
to zero, or or maximal norm, along the whole trajectory. Such kind of trajectories can be 
excluded for practical applications and are thus discarded in the sequel. 

Lemma 3.2. Consider a nondegenerate extremal. Then: 

1. The set T :— {t £ [0,tf] \ p v (t) = 0} has a finite cardinal. 
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2. There exists a measurable function a on [0,tf], with values in [0, 1], such that 

«(0 a.e. on[0,t f }. (12) 

\\Pv{t)\\ 



3. Set *(t) := ^||p»(t)|| -bp m (t). Then 

«(*) = 



i/tf(t)<0, 

1 i/#(t)>0. 

Proof. If f E T, then by the costate equation (fTTj) . p„(t) = —p r (t) is not zero (since the 
extremal is not degenerate). Therefore T has only isolated points, and hence, has a finite 
cardinal. 

Writing «; = ad, with a = ||tu|| and d of unit norm, we get 3>t(iu) = a f^^j (Pv(t), d) — bp m (t) 

Since p v (t) =fi a.e., points 2 and 3 of the lemma follow immediately from the maximization 
condition. □ 



The continuous function W defined in Lemma 13.21 is called switching function. In the 
conditions of the lemma, the extremal control is either equal to 0, or saturating the constraint 
and of direction p v (t). The remaining case, not treated in this lemma and analyzed next, is 
the case where the function vanishes on a (closed) subset I C [0, f/] of positive measure. 

Remark 3.1. Let [to,ti] be a subinterval of / on which a(t) > 0. Then, the control wirt^^] 
is singular. 

Indeed, it suffices to notice that, using ITT 



dw" w/ V^M ; " v 7 lb«(*)ll \\Pv(tW 

and to use the Hamiltonian characterization of singular controls recalled in Section 12.21 

Singular arcs may thus occur in our problem whenever \& vanishes, and we next provide 
an analysis of that case, and show how to derive an expression of such singular controls. 

3.2 Analysis of singular arcs 

Throughout this section, we assume that 

m = -£rAPv(t)\\-bp m (t)=V (13) 

for every t E I, where I is a (closed) measurable subset of [0,t/] of positive Lebesgue 
measure. 

Usually, singular controls are computed by derivating this relation with respect to t, until 
u appears explicitly. The following result is required (see [HI Lemma p. 177]). 
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Lemma 3.3. Let a, b be real numbers such that a < b, and f : [a,b] — > M be an absolutely 
continuous function. Let J be a subset of {t £ [a, 6] | f(t) = 0} os positive Lebesgue measure. 
Then f'(t) — a.e. on J. 

Using this lemma, and extremal equations (fTT|) . one gets, for a.e. t € I, 

bC 

m = ^nW (l|A,( * )III|u(t)l1 ~ <Pv(t)Mt)))+Z(r(t),v(t)Mt),p(t)) -0, (14) 
where the function 

Db . , C ( . (p v ,v) /8D \ 

^{r,v,m,p) = ,, {Pv,v) H r — rr I (p v ,Pr) H n — TT \ "a - '^ / 

m^||u|| m l|P«ll V m \\ v \\ \ I 

dD\\ Vv f D ( Pv ,v) 2 \ 

dmm\\v\\ m \\v\\ 3 J 

does not depend on u. From Lemma f3.2l the relation (fT^|) holds almost everywhere, and 
hence the first term of (|14p vanishes. Therefore, 

*(t)=S(r(t),«(t),m(t),j>(*))=0, (15) 

for almost every < 6 J (actually over every subinterval of positive measure, since the above 
expression is continuous). 

Relations (fT3"l) and (fT4"|) are two constraint equations, necessary for the existence of a 
singular arc. Derivating once more, using Lemma l3.3[ leads to 

*(f) = 0, a.e. on I. (16) 

The control u is expected to appear explicitly in this latter relation. However, since cal- 
culations are too lengthy to be reported here, we next explain how (|16[) permits to derive 
an expression for a(t), and hence, from p2[) . an expression for u{t). When derivating (|15[) . 
the terms where the control u appears are the terms containing v, p m , and m. Recall that 
m = —b\\u\\, that p m — ^{Pv, — ^m"^ ]rf|T + C^)i an d that v is affine in u. Hence, since 
> 0, it is not difficult to see that this derivation leads to an equation of the form 

A(r,v,m,p r ,p v ,p m )a = B(r,v,m,p r ,p v ,p m ), (17) 

almost everywhere on I. This relation should be " generically" nontrivial, that is, the coeffi- 
cient A should not be equal to zero. This fact proves to hold true on numerical simulations. 
We explain below rigorously why this is true generically at least in the case of a scalar 
control (recall that we deal here with a three-dimensional control) . For a scalar control, the 
control system ((l) is of the form 

q = f (q) + uf 1 (q), (18) 
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where /o and f\ are smooth vector fields, and q is the state. In this case, it is well known 
(see e.g. [5]) that, if u is a singular control on /, then there must exist an adjoint vector p 
such that 



The situation encountered here for 3D Goddard's problem is similar to that case: Equa- 
tions (HU), (HO]), (|2T]). are respectively similar to Equations (flB"]). (fT5]) . (fTT)]) ; Equations (jT^|) . 
(|2TJ|) (similarly, Equations (fT3"|) . (fT5|l ) are constraint equations, and Equation (f^Tj) (simi- 
larly, Equation (116p ) permits in general to derive an expression for the control u. The 
vocable "generic" employed above can now be made more precise: it is proved in [6] that 
there exists an open and dense (in the sens of Whitney) subset Q of the set of couples of 
smooth vector fields such that, for every control system (|18p with (Jo, /i) £ G, the set where 
(p? [/oj vanishes has measure zero, and hence Equation (|2ip always permits to de- 

rive u. Additionaly, we can notice that the classical one-dimensional Goddard problem can 
be formulated as a particular case of the general 3D problem described here. In this case, 
it is well known that the second derivative of the switching function provides the expression 
of the singular control, so we can safely assume that [T7] is nontrivial for the restriction to 
the ID problem. Based on these arguments, we should expect the coefficient A of Equa- 
tion p7|) to be non zero in general. This is indeed the case in our numerical simulations 
presented next. Of course, once a(t) has been determined, one has to check (numerically) 
that < a(t) < 1, so that the constraint ||u|| < 1 is indeed satisfied. Here also, numerical 
simulations show the existence and admissibility of such singular arcs (see Section 0|) . 

3.3 Conclusion 

We sum up the previous results in the following theorem. 

Theorem 1. Consider the optimal control problem of maximizing a final cost g°(tf,x(tf)), 
for the control system (QP, with initial conditions f3J) and final conditions x(tf) £ M\. We 
assume that Assumption (H) holds. Letu be an optimal control defined on [0,i/], associated 
to the trajectory (r(-), i>(-), m(-)). Then, there exist absolutely continuous mappings p r (') '■ 
[0,tf] — ► M 3 , Pv(') '■ [0,tf] — * JR 3 , Pm(') '■ [0,tf] — ► 1R, and a real number p° < 0, such that 
(Px('),Pv('),Pm('),P°) is nontrivial, and such that Equations f77]) hold a.e. on [0,tf]. Define 
the switching function "J on [0,tf] by 



(P, /!(<?)) 

<p, [/o,A(g)]> 

(P, [foAhJl(q)})+u(p, [/i,[/ ,/i(g)]> 



on I, 
on /, 
a.e. on /. 



(19) 
(20) 
(21) 



*(*) 



C 



\\p v (t)\\ - bp m (t). 



m(t) 



Then, 



• iff(t) < then u{t) = 0; 
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. if*(t)>Qthenu{t) = ^^; 

• ifty(t) — on a subset I C [0, tf] of positive Lebesgue measure, then Equation {7 
must hold on I , and 

u(t) = a(t)j — -—tt a.e. on I, 

\\Pv(t)\\ 

where a(t) e [0, 1] is determined by |_?7[ ). 

Remark 3.2. The optimal control is piecewise either equal to zero, or saturating the con- 
straint with the direction of p v (t), or is singular. Notice that, in all cases, it is collinear to 
p v (t), with the same direction. 

Remark 3.3 (Optimality status). The maximum principle is a necessary condition for op- 
timality. Second-order sufficient conditions are usually characterized in terms of conjugate 
points (see e.g. [HE]. Unfortunately standard theories do not apply here for two reasons: 
first, the equation in m(t) involves the term ||u(£)|| which is not smooth; second, the struc- 
ture of trajectories stated in the theorem involves both bang arcs and singular arcs, and up 
to now a theory of conjugate points that would treat this kind of trajectory. 

We mention however below a trick, specific to the form of our system, which permits to 
apply the standard theory of conjugate points on every subinterval J of [0, tf] on which u is 
singular and < ||u(t)|| < 1. Let J be such a subinterval. Then, m/0 a.e. on J, and the 
system can be reparametrized by —m(t). Then, denoting q — (r,v), system (jTJ yields 



Now, set 



dq 1 , ui . u 2 \, u 3, \ 

dm u ff u \u\ 



1 a Ul a a U2 a ■ a U3 ■ a 

v = - — - , and - — - = cos Ui cos O2 , t — n — cos v\ sm O2 , r, — rr = sm "2 , 

\\ u \\ \\u\\ ~ ||«|| \u\ 



and consider as new control the control u = (v, ^1,^2)- Notice that the controls Q\ and 
62 are unconstrained, and that v must satisfy the constraint v > 1. However, along the 
interval J it is assumed that < |[tt(i)[| < 1, and thus v does not saturate the constraint. 
Hence, the standard theory of conjugate points applies and the local optimality status of 
the trajectory between its extremities on J can be numerically checked, for instance using 
the code COTCOT (Conditions of Order Two and COnjugate times), available on the 
web3, developed in 0]. This reference provides algorithms to compute the first conjugate 
time (where the trajectory ceases to be optimal) along a smooth extremal curve, based 
on theoretical developments of geometric optimal control using second order optimality 
conditions. The computations are related to a test of positivity of the intrinsic second order 
derivative or a test of singularity of the extremal flow. 

It can be checked as well that every smooth sub-arc of the trajectory is locally optimal 
between its extremities. However, the problem of proving that the whole trajectory (i.e., a 



1 http:/ /www. n7.fr/apo/cotcot 
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succession of bang and singular arcs) is locally optimal is open. Up to now no conjugate 
point theory exists to handle that type of problem. Of course, one could make vary the times 
of switchings but this only permits to compare the trajectory with other trajectories having 
exactly the same structure. A sensitivity analysis is actually required to treat trajectories 
involving singular subarcs. 

4 Numerical experiments 

In this section, we provide numerical simulations showing the relevance of singular arcs in 
the complete Goddard's Problem. For given boundary conditions, the optimal trajectory 
is first computed using indirect methods (shooting algorithm) combined with an homotopic 
approach. Then we use a direct method (based on the discretization of the problem) to 
check the obtained solution. All numerical experiments were led on a standard computer 
(Pentium 4, 2.6 GHz). 

4.1 Numerical values of the parameters of the model 

We implement the optimal control problem of maximizing m(tf) for the system (JTJ) , with 
the constraint ^j. The equations of motion can be normalized with respect to r(0), m(0), 
and go- We follow [15) (in which 2D-trajectories with maximization of the final velocity are 
studied), and set the following parameters. 

• The distance unit is the Earth radius Rt = 6378 10 3 m. 

• Maximal thrust modulus C = 3.5; 6=7. 

• Gravity g(r) = ijfp-r, with g = 1. 



4.2 Numerical simulations with indirect methods 

In our simulations presented hereafter, we prefer to express the objective of the optimal 

control problem in the following form. 

Maximizing m(tf) is equivalent to minimizing the cost 



• Drag D(r,v) = if D |H| 2 e- 500 (H r ll- 1 ) with K D = 310. 



• Initial and final conditions 



r = (0.999949994 0.0001 0.01), v = (0 0), to = 1, 
Tj = (1.01 0), Vf is free, mj is free. 
tf is free. 
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and we assume that there are no minimizing abnormal extremals, therefore the adjoint vec- 
tor can be normalized so that p° = — 1. The results of the simulations are consistent with 
this assumption. 

According to Section 12.11 the Hamiltonian of the optimal control problem under consid- 
eration is 

H = (Pr, V) + ( Pv , ^R^AJL _ g{r) + C— \ - (1 + bp m )\\u\\, 

\ m \\v\\ m/ 

The only difference with the Hamiltonian in l2.1l for the Max m(tf) objective is the additional 
"—1" in the ||u|| term, which leads to the switching function ijj{t) — -^^\\p v (t)\\~ (l+bp m (t)) , 

• if ijj(t) < then u(t) = 0; 

. if > then = ^J^; 

• if ip(t) — on / C [0, tf], then Equation (fl5|) must hold on /, the control u is singular, 
and 

u(t) — a(t) J )v \ \r. a.e. on /, 

\\Pv{t)\\ 

where a(t) G [0, 1] is determined by (fTT|) . We check numerically that < a(t) < 1. 

Furthermore, on a singular subarc, derivating the switching function twice yields the 
expression of a via a relation of the form A(x 1 p)a = B(x,p), see 1 171 The computations are 
actually quite tedious to do by hand, and we used the symbolic calculus tool Maple. The 
expressions of A and B are quite complicated and are not reported here. 

The free final time problem is formulated as a fixed final time one via the usual time 
transformation t = tf s, with s £ [0, 1] and tf an additional component of the state vector, 
such that tf — and t/(0),i/(l) are free, with the associated costate satisfying p\, = —H. 

Transversality conditions on the adjoint vector yield p v (l) = (0 0), p m {X) = 0, and 
pt / (0)=jx / (l)=0. 



4.2.1 Homotopic approach 

In the indirect approach, it is necessary to get some information on the structure of the 
solutions, namely, to know a priori the number and approximate location of singular arcs. 
To this aim, we perform a continuation (or homotopic) approach, and regularize the original 
problem by adding a quadratic (||m|| 2 ) term to the objective, as done for instance in [T31ll9j. 
The general meaning of continuation is to solve a difficult problem by starting from the 
known solution of a somewhat related, but easier problem. By related we mean here that 
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there must exist a certain application h, called a homotopy, connecting both problems. Here, 
we regularize the cost function by considering an homotopic connection with an energy, 



where the parameter of the homotopy is A € [0, 1]. The resulting perturbed problem (P\) 
has a strongly convex Hamiltonian (with respect to u), with a continuous optimal control, 
and is much easier to solve than (P) — (Pi). Assuming we have found a solution of (Po), 
we want to follow the zero path of the homotopy h until A = 1, in order to obtain an 
approximate solution of (P) (or at least sufficient information). The continuation can be 
conducted manually, by finding a suitable sequence (A&) from to 1. However, finding such 
a sequence can be quite difficult in practice, which is why we chose here to perform a full 
path-following continuation. Extensive documentation about path following methods can 
be found in [S]. We use here a piecewise- linear (or simplicial) method, whose principle is 
recalled briefly below. The reason behind the choice of this method over a more classical 
predictor-corrector continuation (such as detailed for instance in [9]) is that we expect the 
problem to be ill-conditioned, due to the presence of singular arcs, which is indeed the case 
in the numerical experiments. 

Simplicial methods. PL continuation methods actually follow the zero path of the ho- 
motopy h : R" +1 — > R n by building a piecewise linear approximation of h. The search 
space R n+1 is subdivided into cells, most often in a particular way called triangulation in 
simplices. This is why PL continuation methods are often referred to as simplicial methods. 
The main advantage of this approach is that it imposes extremely low requirements on the 
homotopy h: since no derivatives are used, continuity is in particular sufficient, and should 
not even be necessary in all cases. 

Definition 4.1 (Simplices and faces). A simplex is the convex hull of n + 1 aflmely inde- 
pendent points (called the vertices) in R n , while a fc-face of a simplex is the convex hull of 
k vertices of the simplex (k is typically omitted for n-faces, which are just called faces, or 
facets) . 

Definition 4.2 (Triangulation). A triangulation is a countable family T of simplices of R™ 
such that the intersection of two simplices of T is either a face or empty, and such that T is 
locally finite (a compact subset of R n meets finitely many simplices). 

Definition 4.3 (Labeling). We call labeling a map I that associates a value to the vertices 
Vi of a simplex. We label here the simplices by the homotopy h: l(v % ) — h(z l ,X l ), where 
v l = (z 1 , A 1 ). Affine interpolation on the vertices thus gives a PL approximation hx of h. 

Definition 4.4 (Completely labeled face). A face [v±, ..,V n ] of a simplex is said completely 
labeled if and only if it contains a solution v e of the equation hx(v) = e = (e, .., e n ), for every 
e > sufficiently small. 




(22) 
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Figure 1: Illustration of some well known triangulations of R x [0, 1] ([0, 1[ for J3): Frcudcn- 
thal's uniform K\, and Todd's refining J4 and J3 



Lemma 4.1 (0 Chapter 12.4]). Each simplex possesses either zero or exactly two completely 
labeled faces (called a transverse simplex in the latter case). 

The constructive proof of this property, which gives the other completely labeled face of 
a simplex that already has a known one, is often referred to as PL step, linear programming 
step, or lexicographic minimization. Then there exists a unique transverse simplex sharing 
this second completely labeled face, that can be determined via the pivoting rules of the 
triangulation. 

A simplicial algorithm thus basically follows a sequence of transverse simplices, from a 
given first transverse simplex with a completely labeled face at A — 0, to a final simplex 
with a completely labeled face at A = 1 (or 1 — e for refining triangulations that never reach 
1), which contains an approximate solution of h(z, 1) = 0. 

4.2.2 Preliminary continuation on the atmosphere density 

In our case, even solving the regularized problem (Pq) is not obvious, due to the aerodynamic 
forces (drag). For this reason, we introduce a preliminary continuation on the atmosphere 
density, starting from a problem without atmosphere. Technically, this is done by using an 
homotopy with the modified parameter 

K e D =9K D , 6 [0,1], 

where Kjj appears in the model of the drag. The shooting method for the problem without 
atmosphere at 8 = converges immediately with the trivial starting point zq — (0.1 0.1 0.1 0.1 0.1 0.1 0.1). 
We would like to emphasize the fact that we have here no difficulties to find a starting point 
for the shooting method. The path following is then achieved with an extremely rough 
integration formula (Euler with only 25 steps), since we just seek a starting point for the 
main homotopy. Thanks to the robustness of the simplicial method, we can afford such a 
low precision to save computational time. The border at 6 = 1 is reached after crossing 
about 120 000 simplices, for a CPU time of 48 seconds. 

Remark 4.1. The adaptive meshsize algorithm described in |14] here strongly reduces the 
oscillations along the zero path, as shown on Figure [2j which decreases the number of 
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simplices required to reach 9 = 1. We can see that the path following using a fixed uniform 
meshsize actually converges to another point, which corresponds to an incorrect solution 
(the final condition on T2 is not satisfied). 




Z5 Z6 Z7 ZS 



Figure 2: Path following for the atmosphere homotopy: fixed uniform triangulation (grey) 
and adaptive meshsize (black). 

The solution we obtain is sufficient to initialize the shooting method at the beginning of 
the main homotopy. Figure [3] represents the solutions of the regularized problem (Pq) for 
9 = and 9 = 1, i.e., without atmosphere and with a normal atmosphere. 

Notice that a direct continuation on the atmosphere with the original non regularized 
problem (P) fails. During the continuation, the process abruptly diverges at a certain value 
for 9, certainly due to the appearance of the singular arc. 

4.2.3 Main continuation on the quadratic regularization 

We now perform the main continuation on the cost (|22p . Figure[3]represents the solutions for 
A = 0, 0.5 and 0.8. It is visible that this continuation process permits to detect the singular 
structure of the solution. The shape on the switching function and of the control norm 
graphs are particularly interesting concerning suspicion of singular arcs. Indeed, we observe 
that, on a certain time interval (roughly [0.1,0.4]), the switching function comes closer to 
zero as A increases, while the control norm keeps values in (0,1). Along the solution for 
A = 0.8, we can guess the appearance of a small arc where ||u|| = 1 at the beginning. These 
facts strongly suggest the appearance of a singular arc. 

With a fixed meshsize of 10~ 4 , the path following takes about 900000 simplices and 
350 seconds to reach A = 0.8, again with an extremely rough integration (Euler, 25 steps). 
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Figure 3: Regularized problem (Po): solutions for 9 — (no atmosphere) and (9=1 (normal 
atmosphere). 

Trying to go further becomes extremely difficult since we lose the singular structure and 
encounter trajectories with incorrect bang-bang structures. However the knowledge of the 
solution for A = 0.8 happens to be sufficient to solve the problem: it provides a good starting 
point for which the shooting method applied to the original problem (P) converges. 

Remark 4.2. This path is more difficult to follow than the previous one for the atmosphere 
homotopy, and the adaptive meshsize algorithm does not work well. We thus use a fixed 
meshsize to perform this homotopy. 

4.2.4 Shooting method applied to the original problem (P) 

When implementing a shooting method (see for instance [21 [El EH HI] ) ; the structure of the 
trajectory has to be known a priori. The structure of the control must be prescribed here by 
assigning a fixed number of interior switching times that correspond to junctions between 
nonsingular and singular arcs. These times {ti)i = i..n smitch are part of the shooting unknowns 
and must satisfy some switching conditions. Each arc is integrated separately, and matching 
conditions must be verified at the switching times, as drawn on the diagram below. 

Unknown: z 

| IVP unknown at to \ (x L ,p L ) \ ... \ (x s ,p s ) \ t± \ ... \ i~ a | 

Value: S Si ng(z) 
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Figure 4: Main homotopy - solutions for A = 0, 0.5 and 0.8. 



f Switch cond (t{) | Match cond (t- L ) \ ... \ Switch cond (t s ) \ Match cond {t s ) | TC(fj) | 

Here, matching conditions reduce to imposing state and costate continuity at the switch- 
ing times. 

A switching condition indicates a change of structure, which corresponds here to an 
extremity of a singular arc. Along such a singular arc, it is required that ip = ip = 0. The 
control is computed using the relation -0 = 0. Therefore, using this expression of the control, 
switching conditions consist in imposing cither ijj = at the extremities of the singular arc, 
or ip — ip — at the beginning of the arc. In our simulations, we choose the latter solution 
which happens to provide better and more stable results. 

The previous results, obtained with an nomotopic approach, provide an indication on 
the expected structure of the optimal trajectory for the original problem (P). Inspection 
of Figure |3] suggests to seek a solution involving a singular arc on an interval [ti,t2]j with 
to < t\ < t2 < tf. As a starting point of the shooting method, we use the solution previously 
obtained with the homotopy on the cost at A = 0.8. 

The IVP integration is performed with the radau5 code (see [H]), with absolute and 
relative tolerances of, respectively, 10 -6 and 10~ 6 . The shooting method converges in 17 
seconds, with a shooting function of norm 5 10 . The condition number for the shooting 
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function is quite high (about 10 12 ), which was expected. The overall execution time of the 
whole approach (preliminary atmosphere homotopy, regularization homotopy, final shoot- 
ing) is about 400 seconds. 

At the solution, the free final is 0.2189, and the objective value is 0.3994, which corre- 
sponds to a final mass of 0.6006. The evolution of altitude, speed and mass during the flight 
are represented on Figure 

ALTITUDE tl|r||— ||r ||) SPEED |||v||) MASS 




Figure 5: Solution with singular arc: altitude, speed and mass. 

We show on Figure [H] the control and switching function. The singular arc is clearly 
visible on the control norm graph. 

4.3 Numerical simulations with direct methods 

In order to validate the solution obtained previously with the shooting algorithm, we next 
implement a direct method. Although direct methods can be very sophisticated (see for 
instance [51 |23j), we here use a very rough formulation, since our aim is just to check if the 
results are consistent with our solution. We discretize the control using piecewise constant 
functions, and the state is integrated on [0,tf] with a basic fixed step Runge-Kutta fourth 
order formula. The values of the control at the discretization nodes, as well as the final 
time tf, thus become the unknowns of a nonlinear constrained optimization problem, the 
constraints being the final conditions for the state. To solve the optimization problem, we 
use the ipopt solver, which implements an interior point algorithm with a filter line-search 
method (see [53J for a complete description). 

With standard options, the algorithm converges after 193 iterations (and 210 seconds) to 
a solution with a final time of 0.2189 and a criterion value of 0.3997. This solution is clearly 
consistent with the results of the shooting method, as shown on Figure [3 which represents 
the norm of the control for the shooting method solution, the direct method solution, and 
a bang-bang reference solution (see below). 
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Figure 6: Solution with singular arc: control and switching function. 
Comparison with a bang-bang solution 

Recall that the usual launch strategy consists in implementing piecewise controls either 
saturating the constraint or equal to zero. To prove the relevance of the use of singular 
controls in the control strategy, we next modify slightly the formulation above in order to 
find a bang-bang solution. Our aim is to demonstrate that taking into account singular arcs 
in the control strategy actually improves (as expected) the optimization criterion. 

We implement a "on-off" structure, with only one switching time t a g. The control is 
chosen so as to satisfy ||u(t)|| = 1 for t < t < t a g, and u(t) = for t s < t < tf. Here, 
the unknowns of the optimization problem are tf, i ff and the direction of the control at the 
discretization nodes before t a g. We obtain a solution with tj = 0.2105, t a s — 0.0580, and 
the value of the criterion is 0.4061, which represents a loss of about 1.6% compared to the 
solution with a singular arc. On this academic example, the gain of the optimal strategy, 
involving a singular arc, over a pure bang-bang strategy, is quite small. This simplified 
problem is a first step in the study of a realistic launcher problem, and permits to illustrate 
the method. 
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Figure 7: Control norm for the shooting and direct method. 
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